A Worm Algorithm for Two-Dimensional Spin Glasses 
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A worm algorithm is proposed for the two-dimensional spin glasses. The method is based on a 
low-temperature expansion of the partition function. The low-temperature configurations of the spin 
glass on square lattice can be viewed as strings connecting pairs of frustrated plaquettes. The worm 
algorithm directly manipulates these strings. It is shown that the worm algorithm is as efficient 
as any other types of cluster or replica-exchange algorithms. The worm algorithm is even more 
efficient if free boundary conditions are used. We obtain accurate low-temperature specific heat 
data consistent with a form c ^ exp( — 2 J/(fc_gT)) , where T is temperature and J is coupling 
constant, for the ±J two-dimensional spin glass. 
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Spin glasses have been studied for many years from 
various points of views 0]. However, the nature of 
low-temperature phases is still not clarified. Much of 
the work on spin glasses relies on computer simulations. 
Monte Carlo simulation has been one of the main tools. 
The Metropolis algorithm proposed more than half a cen- 
tury ago served well for many of the simulational work, 
but for spin glasses, it is hampered by extremely long re- 
laxation times at low temperatures. There have been 
a number of more efficient techniques, noticeably the 
simulated tempering and paralM temp ering I3j|. In 
two-dimensions (2D), cluster algorithms 0,0,0,01 ^^i^* 
which are quite efficient. 

The idea of generating random walks of loops in Monte 
Carlo simulation has a surprisingly longhistory 0| . The 
loop algorithms for quantum systems j9( are ex amp les. 
Recently, several worm algorithms are proposed [lO|, 
for the ferromagnetic Ising and other models. Such al- 
gorithms have the advantage that the Monte Carlo up- 
dates are purely local, while their effects are global. The 
classical worm algorithms based on high-temperature ex- 
pansion variables such as tanh( Jij/fcsT) do not work 
for spin glasses as the weights would be negative for 
the antifcrromagnetic interactions. In this paper, we 
propose a worm algorithm for the two-dimensional spin 
glasses. Our starting point can be thought of as a low- 
temperature expansion for the partition function. Our 
algorithm turns out as efficient as the best cluster al- 
gorithms 't'I , replica Monte Carlo or replica exchange 
algorithms 3]. Moreover, large systems can be simulated 
since we may only simulate one system at a time. 

In the following, we outline a base algorithm and then 
show how to enhance it by multi-step moves. We discuss 
the efficiency of few variations of the algorithm. We re- 
port extensive simulation results for the low-temperature 
specific heat with periodic and free boundary conditions. 

The weight of a spin-glass configuration is proportional 
to the Boltzmann factor exp(J2{ij) ^ij'^i'^j / (^bT)) , 
where ai = ±1, and the site i is on an L x L square 



lattice with periodic boundary conditions. The summa- 
tion is over the nearest neighbor pairs. For simplicity, we 



consider the ± J spin glass where Ji^ 



-J or — J with 



equal probability, although the method is not limited to 
this model. By multiplying the weight by a configuration- 
independent constant, we can rewrite it in an equivalent 
form: 
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{ij) 
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where, w = exp(— 2_fi'), K — J/{kBT); the variable 
bij = ^(1 — JijO'iO'j/J) represents presence (1) or ab- 
sence (0) of an unsatisfied bond. The bonds live on the 
dual square lattice. Note that the variables 6y are not 
independent. They should be set up in such a way that 
an even number of bonds incident on an unfrustrated pla- 
quette, while an odd number of bonds incident on a frus- 
trated plaquette. A groundstate is one such that all the 
frustrated plaquettes are paired and connected by strings 
with minimum total length [l^. At excited states, closed 
loops of strings can form. 

The worm algorithm directly manipulates these 
strings. The weight, Eq. can be sampled with a 
"worm" if we extend the phase space to include a path of 
a worm, with a moving head at location i, and a tail at a 
fixed location zq. The weight is exactly the same as be- 
fore, except that the parity requirement for the head and 
tail is reversed. That is, a frustrated plaquette at head or 
tail requires an even number of bonds, and unfrustrated 
plaquette an odd number of bonds. The movement of the 
worm must preserve the constraint on the bonds. A valid 
configuration of the original problem is formed when the 
worm traces out a closed loop. 

The sites and bonds in the following refer to the dual 
lattice. The base algorithm of the worm movement with 
a periodic boundary condition is as follows: 

1. Pick a site io at random as the starting point. Set 
i ^ if). 

2. Pick a nearest neighbor j with equal probability, and 
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move it there with probabihty w^~^^^ . If it is ac- 
cepted, flip the bond variable 6^ (0 ^ 1) along the 
way, update i ^ j. 
3. If i is at the same site as and winding numbers 
are even, one Monte Carlo loop is finished (exit and 
take statistics), else go to step 2. 
We define the winding numbers as the algebraic sum of 
displacements ^ <5r divided by the linear size L when 
the head and tail meet. The requirement that the wind- 
ing numbers must be even is due to the constraints of 
the bonds on the dual lattice and the spins on the origi- 
nal lattice. A one-to-two mapping to spin configurations 
is possible and valid from a bond configuration only if 
the worm winds the system an even number of times in 
both directions. However, for free boundary conditions, 
where the spins at the boundary have fewer neighbors, 
the winding number constraints are not needed. Because 
both the bond variables and spin variables are available, 
any desired thermodynamic quantities, such as spin-glass 
susceptibility, can be obtained. The algorithm is nothing 
but a Metropolis sampling on the extended phase space. 
It is ergodic, any state can be reached with nonzero prob- 
ability, except at zero temperature. We also note that 
the difference between a ferromagnetic Ising model and a 
spin glass is only at the initial configuration. In the case 
of the 2D ferromagnetic Ising model, because of dual- 
ity, the same bond configuration can also be interpreted 
as a high-temperature expansion loop configuration at its 
dual temperature. Thus, t he p resent algorithm is exactly 
a dual algorithm of Ref. for the 2D ferromagnetic 
Ising model. 

The basic algorithm can be systematically enhanced 
by applying the N-fold way or "absorbing Markov 
chain" method Q|. Note that we are only interested 
in the configuration when the worm forms a closed path, 
and do not care about the (Monte Carlo) dynamics while 
the worm is making its way to meet the tail. We can 
apply an n-step acceleration if it is n-step away from the 
tail iQ. 

Consider a set A of states in state space, consisting of 
the current state and a collection of neighborhood states 
of the current state. For example, in our application, we 
can consider all the states reachable from current state 
in n — 1 steps of moves or less. We calculate the escape 
probability, P(i^|/i), of exiting A to state v given that it 
is in state /i. Let W^J^ be the one-step transition matrix 
elements of the Markov chain for states within the set A 
from /i to ly, and W^J^ for one-step transition probability 
with /i € A, but V ^ A, where A is the complement of A. 
Then the escape probability is given by 

oo 
k=0 

(2) 

where / is an identity matrix. The total escape proba- 
bility is one, -^'(^Im) = 1- Given the current state 



/i, the state v G A is sampled with the probability 
P(i^|/i). The escape probability can be calculated ex- 
plicitly. Let the current state be called 0, and its 4 
neighbor states by one step move 1 to 4. A new state 
is uniquely specified by a path of the moves defined by 
the base algorithm, given the current state. One-step es- 
cape probability is P(0 — > ^) = doWo^, = 1,2,3,4. 
da — 1/(1 — Woo) is fixed by normalization. This is 
just the original N-fold way of Bortz ct al. 13] . For 
a two-step move, from the current site to an inter- 
mediate site a and reaching ly, the escape probability is 
P(0 I') ^ doWoaWaul (1 - Waa) where do is again 

fixed by normalization. Probabilities for three or more 
steps are slightly more complicated, but can be worked 
out. In this work, we consider n — Q (no acceleration), 
and n = 1 to 4 step accelerations. 

At very low temperatures, it can take an exceedingly 
long time to generate one loop. In this case, it is actually 
correct to interrupt the simulation by setting a fixed up- 
per limit to the number of steps used for each loop. Those 
attempts that exceed the upper limit will be treated as 
rejected moves. 

We measure the performance of the algorithm by its 
correlation times. The correlation times are defined 
through the correlation functions of the overlapping spin- 
glass order parameter: 

^ \ {Q{t + t')Q{t'))-{Q{t'm 
[ {Q{t'Y)-m')? J/ 

where the angular brackets denote average of Monte 
Carlo loop moves t' and the outer square brackets mean 
average over the quenched random couplings Jy . The 
quantity Q is the overlap of the spins of two independent 
configurations, 

Q^IE^'^'I- (4) 

i 

Figure ^ demonstrates the efficiency of the algorithm 
for various sizes and inverse dimensionless temperature 
K — J/{kBT). The results are from linear fits of the 
form In /(i) = —t/r + cma window [r, 3r] . The correla- 
tion functions are very close to a pure exponential with 
c w 0. The top part. Fig. ^a), shows correlation time 
T in units of loop moves, the bottom part (b) shows the 
central processor unit (CPU) times. It is useful to sep- 
arate the effect of the intrinsic dynamics defined by the 
base algorithm from the speed-up due to differences in de- 
tailed implementation. The correlation time in Fig. Q^a) 
is independent which of the N-fold-way is used since the 
N-fold-way preserves the dynamics. They all give the 
same correlation time in units of number of loops gener- 
ated. The bottom part Fig. GJb) gives the actual CPU 
time, to, in microsecond for one loop generation divided 
by the number of spins. This number is only slightly de- 
pendent on system size L. The overall efficiency should 




FIG. 1: (a) Exponential relaxation times in units of loop 
trials of the worm algorithm for system sizes L x L where 
L = 8, 32, 128, and 512, vs. dimensionless inverse temperature 
K = J/(kBT). For comparison, the result of single-spin-flip 
Metropolis algorithm on L = 128 is also plotted, (b) CPU 
time (/isec) per loop trial per lattice site for a 32 x 32 lattice 
on a 2.4 GHz AMD Opteron computer, without (0) or with 
n-step acceleration. 



be measured by the product of the two. It is interesting 
to note that as K increases, the correlation time r sat- 
urates and becomes K independent. Unfortunately, the 
time it takes for generating one loop increases exponen- 
tially. Comparing different versions of N-fold-way accel- 
eration, we found that there is a big improvement going 
from the base algorithm to one-step N-fold-way. Further 
increasing the step size does not lead to big improvement 
until very low temperatures. It is clear that, at T = 0, 
one step or two step N-fold-way will not be ergodic, the 
system can be trapped in a configuration. However, if we 
allow for sufficiently long-ranged multi-step attempts, we 
can still make moves even at T = 0. 

What is more remarkable is that a free boundary condi- 
tion will make the algorithm even more efficient. This is 
because we can always start from the boundary, or "out- 
side" the system, and the worm will cut the system into 
pieces. The worm has a much easier job terminating as 
there are order L more possibilities to hit the boundary 
comparing to a single site. Conceptually, we can view 
the system as a planar graph. Various ways of exiting 
the system can be regarded as moving to the single dual 
site representing the outside. Fig. |21 gives the correla- 




FIG. 2: Correlation times, rto, measured in actual CPU 
times for single-spin-flip, 1-step N-fold-way worm algorithm 
with periodic boundary condition and with a free boundary 
condition, for L = 128. The slopes a of the straight line fits 
to exp{aK) axe 14.0, 6.5, and 6.2, respectively. 



tion times measured in real CPU time, the quantity tIq. 
The meaning of this quantity is the amount of CPU time 
needed in order to decorrelate the system such that the 
correlation function is reduced to about e~^. This is a 
fair comparison between completely different algorithms, 
but the results depend on the detail implementation of 
the computer programs. From this plot, we see that the 
free boundary condition case is about 20 times more effi- 
cient than periodic boundary conditions. The correlation 
^of the 2D ±J spin glass diverges as exp(2/f) 
If we define the dynamical critical exponent z 
as rto C^, we have z ~ 7.0 and 3.2, for the Metropolis 
single spin flip and the worm algorithms, respectively. 

We now turn to the specific heat of the spin glass at 
low temperatures. In 1988, Wang and Swendsen found 
by their replica Monte Carlo algorithm that the specific 
heat of the 2D ± J spin glass approaches zero according 
to c ~ exp(— 2if ) Since the energy gap from the 
ground states to the first excited states is 4J, we might 
expect that the specific heat should go as exp{—4:K). 
Wang and Swendsen gave an argument in analogous to 
the ID Ising model with periodic boundary condition, 
where although the minimum excitation is also 4J, the 
configuration appears in the form of a pair of kinks, each 
one of them can move freely, so only a single kink with 
energy 2J should be considered "elementary." Indeed, 
for ID Ising model, the specific heat goes as exp{—2K) 
in the thermodynamic limit. 

However, the above interpretation has been challenged 
in Ref. 15] with a method of exact calculation of the 
partition function and recently in jl^ by Monte Carlo 
simulation, but supported in 18J. Thus, the problem is 
still controversial. 

The present algorithms are well-suited to simulate spin 
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FIG. 3: 2D ±J spin-glass reduced specific heat c/K^ vs. K. 
The legend specifies the lattice size L (/ for free bound- 
ary conditions), Monte Carlo steps, and number of samples 
of random couplings. The straight lines are fit to c/K"^ « 
beiqy{~2K) or tiexp(— 4/^), with b indicated on the lines. The 
errors at the largest K are mainly due to sample fluctuations. 



glasses at low temperatures, particularly the free bound- 
ary condition version. However, as we can see from Fig. 12 
the worm algorithms still have difficulty equilibrating the 
system when coupling K > Z. Therefore, we have used a 
combination of several different algorithms. Each Monte 
Carlo step consists of several worm loop moves, one sweep 
of Metropolis single spin flip, replica Monte Carlo be- 
tween systems at neighboring temperatures and at the 
same temperature. Three replicas for each temperature 
are used. With these algorithms, we were able to equi- 
librate the system down to « 6 for L = 8 and 32, or 
/s: w 4 for L = 128. 

In Fig. 13 we present the reduced specific heat cK^'^ 
vs. K . The asymptotic slopes for large K should resolve 
the issue of 2 J vs. 4J controversy. As can be seen from 
the figure, with periodic boundary conditions, the spe- 
cific heat eventually settles to e-xjp{—AK) because of the 
gap of excitations. However, as L increases, the crossover 
to a fast decay appears at a larger and larger value of 
K. According to Ref. 18], the cross-over length scale 
is I ~ exp(if). This "finite-size-effect" length scale ap- 
pears different from the correlation length ^, which goes 
as exp(2A'). On the other hand, the free boundary con- 
dition results have a rather different finite-size depen- 
dence. We see that the asymptotic form of exp(— 2i4r) 
is approached much faster in this case. Of course, since 
the energy gap with free boundary conditions is 2 J, the 
final form must be exp(— 2i4r) for finite sizes. Our data 
indicate that this is also the form in the thermodynamic 
limit. Some weak size dependence on the amplitudes is 
seen, but surprisingly, the thermodynamic limit is ob- 
tained much faster with free boundary conditions than 
with periodic boundary conditions. 

In summary, we have shown numerically that the low- 



temperature specific heat goes like exp(— 2i^) in the 
thermodynamic limit. This is observed much clearer with 
free boundary conditions. Our new calculations support 
the original argument of Wang and Swendsen regarding 
the asymptotic form of heat capacity for 2D ±J spin 
glass. We demonstrate these results with an efficient 
worm algorithm. The worm algorithm presented here 
should be applicable to any models defined on a (planar) 
graph where the concept of dual graph can be defined. 
It does not seem possible for 3D lattices. It is also in- 
teresting to study the clusters generated in the worm 
algorithms and to relate them to other quantities. More 
detailed study will be presented elsewhere. 

This work is supported in part by a Faculty Research 
Grant of National University of Singapore. 
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